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Introduction 

Cavitation  deals  with  the  vaporization  of  a  liquid  as  it  flow  past  or  through  a  body  of 
interest.  As  the  flow  velocity  increases  and  the  local  pressure  decreases,  the  liquid  can 
change  to  vapor  if  the  overall  pressure  is  sufficiently  close  to  the  vapor  pressure.  If  the 
pressure  change  in  a  flowing  liquid  brings  the  liquid  below  its  vapor  pressure,  local 
regions  of  the  fluid  can  change  phase  and  become  vapor.  In  general,  the  phenomenon  of 
cavitation  increases  as  the  overall  pressure  in  the  system  is  reduced. 

While  cavitation  remains  a  very  fascinating  physical  phenomenon,  it  also  is  of  much 
interest  in  engineering  circles  because  of  the  rapid  damage  cavitation  can  do  to  even  the 
toughest  materials.  In  the  presence  of  cavitation,  rapid  pitting  and  other  surface  and 
structural  damage  can  occur  very  quickly.  For  this  reason  improved  understanding  is 
needed  to  circumvent  its  damaging  nature. 

Cavitation  occurs  in  many  different  forms.  In  some  cases,  cavitation  occurs  as  a  bubble 
cloud  moving  through  the  flowfield.  Such  cavitation  is  referred  to  as  cloud  cavitation.  In 
other  cases,  cavitation  appears  as  a  bubble  attached  to  the  surface.  This  phenomenon  is 
known  as  sheet  cavitation.  A  third  type  of  cavitation  is  vortex  cavitation  in  which 
cavitation  bubbles  appear  in  the  low-pressure  core  of  a  propeller  tip  vortex.  The  present 
report  focuses  on  the  second  of  these  types,  attached  or  sheet  cavitation. 

The  dynamics  of  attached  sheet  cavitation  are  highly  complex.  The  flowfield  is  two 
phase  with  both  liquid  and  vapor  being  present.  In  nearly  every  case,  the  flowfield 
contains  significant  unsteady  effects  and  the  flow  is  generally  three-dimensional.  The 
multi-phase  nature  and  the  unsteadiness  both  increase  the  difficulty  of  making  diagnostic 
measurements.  The  presence  of  the  vapor  phase  can  make  the  fluid  opaque  so  that 
optical  diagnostics  cannot  be  employed.  In  addition,  the  unsteadiness  makes  it  difficult 
to  document  the  local  nature  of  the  flow.  Despite  these  difficulties,  large  numbers  of 
experiments  have  been  conducted,  and  there  is  a  considerable  volume  of  data  available  in 
the  literature.  Certainly,  this  existing  data  has  given  much  light  on  cavitation  processes, 
but  detailed  local  understanding  of  the  cavitation  phenomena  is  still  lacking.  The  present 
study  looks  at  applying  detailed  CFD  computations  to  cavitating  flowfields  as  an  aid  for 
guiding  and  improving  experiments  and  as  a  means  of  providing  increased  understanding 
of  the  experimental  results  obtained  to  date. 


Experimental  observations  of  attached  cavitation  have  provided  us  with  a  considerable 
understanding  of  cavitation  processes.  Experimental  observations  of  cavitation  on  a 
hydrofoil  have  consistently  shown  that  the  surface  pressure  lies  below  the  vapor  pressure. 
In  addition,  the  location  at  which  cavitation  appears  generally  exhibits  a  pressure  that  is 
somewhat  below  the  vapor  pressure  for  the  fluid  temperature  of  interest.  The  cavity  can 
appear  in  the  nose-region  of  a  hydrofoil,  at  the  mid-chord  position,  or  near  the  trailing 
edge.  A  particularly  appropriate  study  for  the  present  analysis  is  the  two-part  study  of 
Franc  and  Michel  (1996).  Their  data  clearly  shows  that  the  location  at  which  cavitation 
appears  on  an  isolated  hydrofoil  is  a  function  of  the  angle  of  attack  as  well  as  the 
cavitation  number. 

Attached  cavitation  typically  produces  long,  thin  cavitation  regions.  The  aft-end  of  the 
bubble  is  generally  highly  unsteady  and  is  composed  of  a  two-phase  region  in  which  the 
fluid  in  the  vapor  is  converted  back  to  liquid.  The  interface  between  the  liquid  and  the 
vapor  can  be  shiny  or  diffuse.  A  shiny  interface  generally  implies  that  the  interface  is 
laminar,  while  a  diffuse  interface  is  generally  characteristic  of  the  presence  of  turbulence. 
Of  particular  interest  is  the  experimental  observation  that  the  leading  edge  of  the  cavity 
can  sometimes  appear  downstream  of  the  minimum  pressure  point  on  a  hydrofoil.  This 
implies  that  the  liquid  successfully  negotiates  the  lowest  pressure  region  in  the  flowfield 
without  changing  to  vapor,  but  then  cavitates  (or  switches  to  vapor)  at  some  higher 
pressure.  This  clearly  suggests  that  non-equilibrium  phenomena  can  be  a  part  of  the 
cavitation  process. 

Cavitation  Models 

Over  the  years  a  number  of  cavitation  models  have  been  used  to  study  the  cavitation 
process.  The  simplest  of  these  treats  the  cavity  as  a  constant  pressure  region  that  is  filled 
with  vapor.  The  pressure  in  the  constant  pressure  cavity  region  is  set  equal  to  the  vapor 
pressure  of  the  liquid.  The  interface  between  the  vapor  and  the  liquid  is  then  treated  as  a 
constant  pressure  streamline.  The  shape  of  the  interface  is  dictated  by  the  dynamics  of 
the  liquid  flow  and  the  trajectory  of  the  constant  pressure  interface.  Since  the  interface  is 
a  streamline,  it  is  clear  that  in  this  ‘cavitation’  model,  no  cavitation  occurs.  The  liquid 
flows  around  the  bubble  and  the  pressure  of  the  liquid  dictates  the  shape  of  the  interface. 
The  cavity  consists  of  a  fixed  volume  of  vapor,  but  the  vapor  inside  this  bubble  is 
quiescent.  Consequently,  this  cavitation  model  does  not  involve  cavitation  at  all  (expect 
perhaps  at  very  early  times  when  the  liquid  inside  the  cavitation  bubble  was  initially 
converted  to  vapor). 

A  major  difficulty  with  this  type  of  model  is  that  it  requires  some  ‘closure’  condition  at 
the  aft  end  of  the  cavitation  bubble.  The  closure  condition  can  take  on  several  different 
forms,  but  one  classical  method  is  to  use  an  artificial  body  or  wall  to  guide  the  liquid 
back  from  the  interface  (which  will  normally  be  above  the  hydrofoil  surface)  so  that  it 
becomes  attached  to  the  hydrofoil  surface  once  more.  Combination  of  this  simple 
cavitation  model  with  a  velocity  potential  model  or  a  full  Navier-Stokes  formulation  is 
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generally  quite  effective.  The  constant  pressure  cavitation  region  gives  reasonable 
agreement  with  experiment,  but  it  does  lack  fundamental  understanding. 

The  emphasis  in  the  present  research  is  upon  using  a  more  physically  realistic  description 
of  cavitation  that  can  help  to  provide  more  mechanistic  understanding  of  cavitation 
phenomena  and  can  lead  to  improved  understanding  and  control. 


Present  Approach 

The  modeling  approach  used  in  the  present  research  is  based  upon  highly  resolved 
numerical  solutions  of  the  full  Navier-Stokes  equations  coupled  with  a  detailed  model  for 
the  presence  of  cavitation.  To  provide  an  appropriate  vehicle  for  incorporating  the 
cavitation  model,  the  Navier-Stokes  equations  are  expressed  in  a  generalized  form  that 
describes  an  arbitrary  fluid  with  an  arbitrary  equation  of  state.  The  arbitrary  fluid 
description  facilitates  the  extension  to  a  multi-phase  model  of  cavitation. 

To  provide  appropriate  background,  we  start  by  writing  the  Navier-Stokes  equations  for  a 
pure  (single-phase,  single-component)  fluid.  We  then  introduce  the  auxiliary  cavitation 
model.  Although  most  of  our  results  are  based  upon  a  two-fluid  model,  we  briefly 
summarize  a  single-fluid  model  that  we  have  also  tested.  Both  the  single-fluid  model  and 
the  two-fluid  model  follow  directly  from  the  generalized  Navier-Stokes  formulation. 

The  Navier-Stokes  Equations:  In  vector  notation,  the  equations  of  motion  can  be  written 
as: 


dQP  dQ  dE  dF  dG 

r — —+— + — + — + — 

dr  dt  dx  dy  dz 


=  v.t. 


In  this  expression,  ( x,y,z )  represents  the  Cartesian  coordinates,  t  is  the  physical  time 
describing  unsteady  transients,  and  x  represents  a  pseudo-time  that  is  used  for  iteration  at 
a  given  time  step.  The  quantities  Qp,Q,E,F,andG  are  vectors  containing  the  primitive 

variables,  the  conservative  variables,  and  the  conserved  fluxes  in  the  x,  y,  and  z  directions 
respectively.  They  are  given  as, 


V 

r  P  ' 

f  P «  ] 

pv  ^ 

r  pw  ^ 

u 

p  u 

p  w2  +  p 

.  p  MV 

p  UW 

V 

Q  = 

pv 

E  = 

P  MV 

F  = 

P  v2+p 

G  = 

pvw 

w 

pw 

p  UW 

pvw 

pw2  +  p 

Jj 

{ph°-p, 

O 

•si 

a 

Q. 

,  f>vh°  1 

"O 

>• 

o 

^  ■ 

3 


The  quantity,  T,  that  multiplies  the  pseudo-time  derivative  is  a  matrix  that  is  used  to 
ensure  that  all  terms  in  the  equations  are  properly  ordered  so  that  the  pseudo-time 
iteration  converges  efficiently  under  various  limiting  conditions. 

The  quantity,  v.t. ,  signifies  the  diffusion  or  ‘viscous’  terms  which  are  given  by, 
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The  vector  of  dependent  variables,  Qp ,  that  appears  inside  the  diffusion  terms  is  the 

same  variable  that  is  used  as  the  primary  dependent  variable  on  the  left-hand  side  of  the 
equation.  Note  that  the  vector  that  appears  inside  the  physical  time  derivative  is  Q  rather 

than  Qp .  The  change  from  Q  to  Qp  simplifies  the  computation  without  compromising 

the  global  conservation  advantages  of  the  conservative  flux  terms. 

The  quantities,  R^.R^,  etc.  that  appear  in  the  diffusion  terms  are  matrices  that  contain 

the  diffusion  properties  of  the  fluid  in  question.  In  most  cases  of  interest  for  cavitation 
modeling,  the  primary  diffusion  properties  are  the  viscosity,  p,  and  the  thermal 
conductivity,  k.  The  physical  properties  of  the  fluid  are  the  primary  terms  in  the  subject 
matrices.  In  general,  all  nine  of  these  matrices  are  highly  sparse.  For  a  general 
compressible  fluid  for  which  the  Stokes  approximation  has  been  made,  these  matrices 
take  the  form: 
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Note  that  the  first  row  of  all  nine  of  the  matrices,  R^.R^,  etc.  is  entirely  zero.  This 

corresponds  to  the  well-known  fact  that  the  continuity  equation  contains  no  diffusion 
terms.  Similarly,  the  first  row  of  all  nine  matrices  is  zero,  indicating  that  the  pressure 
does  not  appear  explicitly  in  the  diffusion  terms.  The  remaining  term  in  the  second,  third 
and  fourth  rows  describe  the  viscous  forces  that  appear  in  the  three  components  of  the 
momentum  equations.  Note  that  only  the  viscosity  appears  in  the  second  through  fourth 
rows.  The  energy  equation  is  the  only  row  that  contains  multiple  entries  in  each  matrix. 
These  multiple  entries  contain  the  heat  conduction  term  (the  diagonal  term  multiplied  by 
the  thermal  conductivity),  and  the  velocity  components  multiplied  by  the  viscous 
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coefficient.  The  product  of  velocity  and  viscosity  represents  the  viscous  dissipation  of 
kinetic  energy  into  heat  energy.  These  terms,  taken  together,  describe  the  so-called 
viscous  dissipation  terms. 

The  equations  of  motion  are  completed  by  an  equation  of  state,  a  thermodynamic  relation 
that  defines  the  enthalpy  in  terms  of  two  other  thermodynamic  variables,  and 
mathematical  functions  describing  the  viscosity  and  thermal  conductivity  as  a  function  of 
two  thermodynamic  variables.  Since  the  temperature  and  pressure  both  appear  in  our 
primary  dependent  variable,  Qp ,  we  choose  to  express  all  four  of  these  parameters  as 

functions  of  the  pressure  an  temperature.  Thus,  the  equation  of  state  is  taken  as  an 
arbitrary  function  that  relates  the  density  to  the  temperature  and  pressure, 

P  =  P(P,T ) 

The  enthalpy  is  likewise  given  in  terms  of  the  pressure  and  temperature  as, 
h  =  h(p,T ) 

The  transport  properties  (viscosity  and  thermal  conductivity)  are  likewise  expressed  as 
arbitrary  functions  of  pressure  and  temperature.  For  the  viscosity, 

H  =  li(p,r) 

and  for  the  thermal  conductivity, 
k  =  k(p,T ) 

Finally  the  stagnation  enthalpy  connects  the  enthalpy  and  the  kinetic  energy, 
h?  =  h  +  i  (w2  +  v2  +  w2  J 


The  Navier-Stokes  equations  given  above,  along  with  these  four  auxiliary  relations,  and 
the  definition  of  the  stagnation  enthalpy  gives  a  fully  defined  system  that  describes  fluid 
flow  under  very  generalized  conditions.  If  the  continuity  equation  is  replicated  so  that  it 
applies  to  multiple  species  or  phases,  the  equations  can  likewise  describe  the  flow  of 
multi-component  and/or  multi-species  flows.  One  of  . our  cavitation  models  uses  a  two- 
phase  description.  Similarly  we  could  include  multiple  momentum  equations  and/or 
multiple  energy  equations  to  cover  these  complexities. 

We  have  tested  two  different  cavitation  models  based  on  this  general  formulation.  One 
uses  a  single-phase  representation,  while  the  other  uses  a  two-phase  description.  These 
two  models  are  described  briefly  below.  As  noted  above,  the  two-phase  model  represents 
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our  preferred  model,  but  its  characteristics  are  more  properly  put  in  focus  by  comparing  it 
with  the  single-phase  fluid. 


Single-Phase  Cavitation  Model:  The  single-phase  cavitation  model  was  the  first  model 
we  used  to  simulate  the  liquid/vapor  phase  change  process.  In  this  model,  the  density  of 
the  fluid  is  treated  as  a  continuous,  single-valued  function  of  the  pressure.  To  account  for 
cavitation,  we  replace  the  discontinuous  change  in  density  at  the  cavitation  pressure  by  a 
rapid,  but  continuous  variation.  For  convenience,  we  also  simplify  the  equation  of  state 
so  that  the  density  is  a  function  of  pressure  only,  p  =  p(p),  instead  of  a  function  of  both 
pressure  and  temperature,  p  =  p (p,T ).  This  allows  us  to  replace  the  energy  equation  by 
the  statement  that  the  temperature  is  constant.  This  approximation  is  not  necessary,  but 
simplifies  the  analysis  slightly  in  that  the  energy  equation  is  uncoupled  from  the 
momentum  equations  for  this  ‘compressible’  fluid  and  need  not  be  solved. 

The  net  result  of  this  single-phase  cavitation  model  is  that  as  the  pressure  is  decreased 
below  the  vapor  pressure  the  density  of  the  water  decreases  from  a  value  corresponding 
to  that  of  a  liquid  to  one  corresponding  to  that  of  a  vapor.  A  transition  in  the  opposite 
direction  occurs  when  the  pressure  is  increased  above  the  cavitation  pressure.  With  this 
model,  the  density  of  the  fluid  in  the  vicinity  of  a  hydrofoil  rapidly  decreases  from  liquid¬ 
like  conditions  to  vapor-like  conditions  in  a  continuous  manner  as  the  surface  pressure 
approaches  the  cavitation  pressure.  Making  the  density-pressure  curve  continuous  rather 
than  discontinuous  as  it  is  in  under  equilibrium  conditions  simplifies  the  numerical 
procedure  while  retaining  much  of  the  physical  characteristics  of  the  cavitation  process. 

In  the  computations,  the  pressure  increment  across  which  the  transition  between  “liquid” 
and  “vapor”  regimes  was  accomplished  was  treated  as  a  parameter.  The  transition 
between  the  two  densities  was  parameterized  by  choosing  a  pressure  increment,  Ap ,  as  a 
fraction  of  the  dynamic  pressure.  Excellent  results  were  obtained  when  A p  was  set  to 
30-40%  of  the  dynamic  pressure,  but  solutions  could  also  be  obtained  when  it  was 
reduced  to  a  few  percent.  Dropping  the  temperature  from  the  equation  of  state  also 
decouples  the  energy  equation  from  the  equations  of  motion  as  in  incompressible 
formulations  although  it  is  emphasized  that  the  present  formulation  contains  a  finite 
speed  of  sound  (in  the  transition  region),  and  allows  large  changes  in  the  density.  This 
single-phase  model  therefore  represents  a  simple  approximation  to  cavitation  in  a 
constant  temperature  fluid. 

In  representative  computations  with  this  model,  the  pressure  increment  across  which  the 
transition  from  the  'liquid'  to  the  'vapor'  regime  was  accomplished  was  treated  as  a 
parameter.  The  transition  between  the  two  densities  was  parameterized  by  choosing  a 
pressure  increment,  Ap ,  as  a  fraction  of  the  dynamic  pressure  and  connecting  the  'liquid' 
density  with  the  'vapor'  density  by  a  cubic  equation  that  was  first-derivative  continuous  at 
either  end.  Excellent  results  were  obtained  when  Ap  was  set  to  30-40%  of  the  dynamic 
pressure,  but  solutions  could  also  be  obtained  when  it  was  reduced  to  a  few  percent. 
Representative  pressure-density  curves  for  the  single-phase  model  are  given  in  Fig.  1  for 
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a  series  of  values  of  A p .  These  curves  give  an  indication  of  the  shapes  of  the  'equations 
of  state'  used  for  the  single-phase  computations. 

As  was  noted  above,  dropping  the  temperature  from  the  equation  of  state  decouples  the 
energy  equation  from  the  equations  of  motion  as  in  incompressible  formulations. 
Nevertheless  it  is  emphasized  that  the  present  formulation  is  fully  compressible.  It  allows 
large  changes  in  the  density  and  results  in  a  finite  speed  of  sound.  For  a  fluid  with  the 
equation  of  state,  p  =  p(p),  the  speed  of  sound  is  given  by  the  square  root  of  the 
pressure-density  derivative, 

_  Mp 
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For  the  present  single-phase  computations,  the  equation  of  state  was  divided  into  three 
sections,  a  vapor  region,  a  transition  region,  and  a  liquid  region.  In  the  vapor  and  liquid 
regions  the  density  was  taken  as  a  linear  function  of  pressure  to  give  a  quantitatively 
realistic  speed  of  sound.  In  the  transition  region,  the  density  was  expressed  as  a  cubic 
function  of  the  pressure.  The  three  regions  are  given  as, 


p  =  c,p 


n<  „  AP 
P^Pv~  — 


P  =  Pv+(P/-PvJ 
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cv  cl 


Vj-vf,  Ij-  -  ' 
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P-Pv 

{  A p 
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py-—&pip,+— 


P  =  cIp  p>pv+-2- 

In  these  expressions,  pi  and  pv  represent  the  pressure  transition  points,  p;  and  pv  are 
the  corresponding  densities,  and  A p  =  Pi  -  pv  signifies  the  width  of  the  transition  region. 
The  quantities,  c/  and  cv  are  the  speeds  of  sound  in  the  liquid  and  vapor.  These  values 
are  taken  from  the  literature.  The  liquid  speed  of  sound,  c/  was  taken  as  1510  m/s, 
corresponding  to  conditions  in  water  at  293K.  The  speed  of  sound  in  vapor,  cv ,  was 

chosen  as  420  m/s  corresponding  to  the  speed  of  sound  in  a  perfect  gas  with  a  molecular 
weight  of  18,  and  a  specific  heat  ratio  of  1.286  at  293K  and  one  atmosphere. 

The  proportionality  constants  in  the  vapor  and  liquid  regions  in  the  example  above  are 
chosen  to  give  a  constant  speed  of  sound  in  the  pure  media  that  is  representative  of  that  of 
the  physical  fluid.  For  computational  purposes,  this  choice  is  arbitrary.  The  common 
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alternative  of  infinite  speeds  of  sound  (incompressible  fluid)  in  each  phase  is  equally 
valid.  In  the  vapor  region,  the  speed  of  sound  should  vary  with  local  conditions,  but  a 
constant  is  appropriate  for  constant  temperature  vapor. 

The  speed  of  sound  in  the  transition  region  varies  continuously  from  the  value  in  the  pure 
vapor  region  to  the  value  in  the  pure  liquid  region,  but  it  reaches  a  quite  low  minimum  in 
the  intervening  region.  In  particular,  since  the  derivative  of  the  pressure-density  function 
is  continuous  at  the  liquid  and  vapor  'boundaries',  the  speed  of  sound  at  either  end  of  the 
transition  region  is  equal  to  that  in  the  appropriate  pure  fluid.  In  the  transition  region,  the 
speed  of  sound  starts  at  the  value  corresponding  to  that  of  the  vapor,  then  decreases 
rapidly  to  a  very  low  value  until  near  the  end  of  the  transition  region  at  which  point  it 
increases  rapidly  until  it  reaches  the  speed  of  sound  in  the  liquid.  This  continuous 
variation  of  the  speed  of  sound  in  a  mixture  with  a  very  low  sound  speed  over  most  of  the 
mixture  region  is  typical  of  vapor-liquid  mixtures.  In  the  present  case,  the  magnitude  of 
the  sound  speed  depends  upon  the  width  (A p  of  the  transition  region.  When  A p  is 
increased,  the  speed  of  sound  in  the  transition  region  increases.  When  A p  becomes  very 
small  the  speed  of  sound  in  the  transition  region  becomes  very  small.  In  the  limit  as  A p 
goes  to  zero  {dp! dp  -»  °°),  the  speed  of  sound  drops  to  zero.  Consequently,  supersonic 
velocities  are  easily  encountered  in  the  vicinity  of  the  transition  region.  Representative 
curves  showing  the  variation  of  the  speed  of  sound  in  the  transition  region  are  given  on 
Fig.  2  for  a  series  of  Ap's. 

Preliminary  results  based  upon  this  model  were  obtained  for  an  ellipse,  a  NACA  66(mod) 
at  one  and  four  degrees  angle  of  attack,  and  for  flow  over  a  cylinder.  Results  for  the 
ellipse  are  given  on  Fig.  3  which  shows  Mach  number  curves  for  flow  over  an  ellipse  at 
minus  2  degrees  angle  of  attack.  Because  the  angle  of  attack  is  negative,  the  high 
velocity  flow  region  occurs  on  the  under  side  of  the  ellipse.  The  speed  of  sound  in  the 
single-phase  liquid  state  is  1520  mis  (corresponding  to  the  sound  speed  in  water),  so  the 

Mach  number  there  is  on  the  order  of  10-3.  Near  the  hydrofoil,  the  flow  passes  through 
the  density  transition  shown  on  Figs.  1  and  2  where  the  speed  of  sound  decreases  rapidly. 
Consequently,  the  Mach  number  in  the  density  transition  region  increases  from  nearly 
zero  to  above  Mach  5  and  then  back  to  the  incompressible  regime  in  a  very  thin  zone. 
Inside  the  high  Mach  number  zone,  the  fluid  is  pure  vapor,  while  outside  it  is  pure  liquid. 
Thus,  the  Mach  number  contours  graphically  describe  the  shape  and  size  of  an  attached 
cavitation  bubble  on  the  underside  of  the  ellipse  near  the  leading  edge. 

Comparisons  between  the  predictions  from  the  single-phase  model  and  experimental  data 
from  Shen  and  Dimotakis  [1]  for  a  NACA  66  hydrofoil  were  also  made,  but  are  not 
shown  herein.  The  comparison  indicated  that  the  single-phase  model  could  replicate  the 
main  features  of  sheet  cavitation  for  both  leading-edge  and  mid-chord  cavitation.  In 
particular  quantities  such  as  surface  pressure  and  cavity  length  were  predicted  reasonably 
well.  A  series  of  different  values  of  A p  (the  width  of  the  density-pressure  transition) 
were  used  to  ascertain  the  effect  on  cavity  length  and  thickness.  For  cases  where  A p  was 
20  -  40%  of  the  dynamic  pressure,  reasonably  good  convergence  to  a  steady  solution  was 
obtained.  When  the  values  of  A p  were  decreased  to  approximately  1%  of  the  dynamic 
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pressure,  the  cavity  interface  became  much  sharper;  however,  the  numerical  stability  of 
the  algorithm  deteriorated  significantly,  even  with  TVD  shock  capturing.  Simulations 
could  be  run,  but  the  values  of  the  time-step  were  too  small  to  obtained  converged 
solutions.  In  assessing  the  results,  it  appeared  that  the  reason  for  obtaining  steady 
solutions  was  numerical  rather  than  physical  in  nature.  At  larger  values  of  Ap  the 
interface  between  the  liquid  and  the  vapor  was  spread  out  over  a  number  of  grid  points 
and  this  appeared  to  stabilize  an  otherwise  unsteady  flowfield.  This  steadiness 
disappeared  as  the  width  of  the  transition  region  was  decreased.  Further,  for  cases  in 
which  the  transition  region  was  too  wide  (i.e.,  the  value  of  Ap  was  increased),  the 
relatively  small  decrease  in  pressure  did  not  allow  the  density  inside  the  cavity  to  drop 
below  the  vapor  density. 

In  addition  to  producing  steady  cavities,  a  second,  and  more  fundamental,  difficulty  with 
the  single  fluid  model  is  that  it  strictly  enforces  an  equilibrium  assumption  on  the 
cavitation  process.  With  the  single-fluid  model,  it  is  impossible  for  the  cavity  to  start 
downstream  of  the  low  pressure  point  on  the  airfoil  (as  is  observed  experimentally).  The 
single-fluid  model  requires  that  cavitation  start  as  soon  as  the  pressure  reaches  the  vapor 
pressure,  so  that  in  all  cases,  the  leading  edge  of  the  cavitation  bubble  will  be  ahead  of 
the  minimum  pressure  point.  This  effect  is  countered  by  the  dual-phase  model  described 
below. 


Dual-Phase  Model:  The  dual  phase  model  treats  the  cavitation  problem  by  considering 
two  distinct  phases,  a  liquid  phase  and  a  vapor  phase.  Each  phase  is  treated  by  a  separate 
continuity  equation,  but  a  single  momentum  equation  is  used  for  the  mixture.  The 
transition  between  phases  is  modeled  by  a  rate  process  which  appears  in  the  separate 
continuity  equations  as  a  source  and  a  sink  of  each  phase.  In  the  present  model,  the 
effects  of  surface  tension  are  ignored  and  both  liquid  and  vapor  are  allowed  to  co-exist 
over  a  range  of  pressures.  The  use  of  a  rate  process  for  converting  liquid  to  vapor  and 
vice-versa  allows  the  cavitation  inception  point  to  occur  downstream  of  the  minimum 
pressure  point  in  agreement  with  experiment  thereby  giving  it  a  meta-stable-like 
character. 

To  incorporate  the  two  phases  into  the  computational  formulation,  the  equations  of 
motion  are  augmented  by  adding  a  second  continuity  equation  and  by  incorporating 
source  terms  to  account  for  mass  conservation  during  phase  change.  For  treating  the  two 
phases,  we  introduce  the  volume  fraction  of  vapor,  or  the  void  fraction,  av,  which 
represents  the  volume  fraction  of  vapor  present  at  any  point.  The  volume  fraction  of 
liquid  is  then  given  as  one  minus  the  volume  fraction  of  vapor  (az  =l-av)  since  the 

sum  of  the  volume  fraction  of  liquid  and  vapor  must  be  unity.  From  the  volume 
fractions,  the  density  of  the  mixture  is  then  given  by 

P  =  «vpv  +  (l-av)pz 
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where  p;  represents  the  density  of  the  liquid  phase,  and  pv  corresponds  to  the  density  of 

the  vapor  phase.  The  model  for  the  two-phase  mixture  uses  a  different  equation  of  state 
for  each  phase.  The  equations  of  state  are  analogous  to  the  single-phase  equation  of  state 
given  above: 

Pz  =P i(p,T),  (liquid phase),  and 
Pv  =  Pv  (p>t  )>  (vapor  phase). 

The  equations  of  motion  of  the  two-phase  system  can  then  be  written  in  a  vector  form 
analogous  to  that  used  above  for  the  single-phase  model. 

Before  going  further,  we  digress  to  compare  the  volume  fraction  with  the  mass  fraction 
and  the  mole  fraction  that  are  traditionally  used  in  mixtures  of  gases.  First  of  all,  note 
that  the  volume  fraction  has  the  units: 

_  volume  of  vapor  _ 
v  volume  of  mixture  n?mix 

The  density  function  described  here  for  the  vapor  (or  for  the  liquid)  has  standard  units  of 
density: 


_  mass  of  vapor  _  kgv 
v  volume  of  vapor  my 

That  is  to  say,  the  density  of  vapor  (liquid)  is  defined  as  the  mass  of  vapor  (liquid) 
divided  by  the  area  that  the  vapor  (liquid)  occupies.  As  a  result,  the  density  of  the 
mixture,  p  =  avpv  +  (l  -  av)p t  has  the  units: 

p  =  avpv  +  (l-av)p/  =  ^  +--H} 

mmix  mv  mmix  ml 

Clearly,  these  are  the  correct  units. 

When  the  mass  fraction,  Yv  is  used,  the  partial  density  is  defined  in  a  different  manner. 
The  mass  fraction  has  the  units, 


_  kgv  ,  kgi  _  kg mix 

3  3  3 

mmix  mmix  mmix 


mass  of  vapor  _  kgv 
mass  of  mixture  kgmix 
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and  when  multiplied  by  the  density  of  the  mixture  we  have, 


p=yvp+(i-rv)p=-fev 


kg 


mix 


kg. 


+ 


kgl  kgmfx  _  kgv  kgi  _  kg 


mix 


mix  m > 


mix 


kg ; 


mix  m, 


mix 


m, 


mix 


ml 


mix 


m, 


two: 


By  comparing  units,  we  can  relate  the  mass  fraction  to  the  void  fraction  as, 


Lmix 


m: 


mix 


ix  kSmix  rt?mix  rt?mix 


From  the  above  relation,  we  can  express  a  ‘partial’  density, 

0  _p Y  _  massqfvapor  _  £gv 
v  v  volume  of  mixture  m^tx 


The  partial  density  corresponds  to  the  density  the  vapor  (liquid)  would  have  if  it  were 
oblivious  to  the  existence  of  the  liquid  (vapor)  and  occupied  the  entire  volume  assigned 
to  the  fluid.  The  partial  density  is  the  density  that  is  typically  used  in  multi-species 
mixtures  of  gases.  Note  we  use  the  superscript  caret  to  distinguish  this  partial  density 
from  the  vapor  density  that  appears  when  we  use  volume  fractions.  Also  note  that  the 
units  are  different.  The  density  of  the  vapor,  pv ,  is  the  mass  of  vapor  divided  by  the 

volume  the  vapor  occupies.  The  partial  density,  pv ,  is  the  mass  of  vapor  divided  by  the 
total  volume  the  fluid  occupies.  In  the  volume  fraction  expression,  the  density,  pv, 

represents  the  standard  density — the  total  mass  of  vapor  divided  by  the  volume  it 
occupies.  In  the  mass  fraction  representation,  the  partial  density  denotes  the  total  mass  of 
vapor  divided  by  the  entire  volume  occupied  by  the  mixture.  This  density  implies  that 
the  vapor  does  not  ‘see’  the  liquid  and  responds  just  as  though  it  were  the  only  specie 
(phase)  present.  The  difference  between  these  two  densities  is  exactly  the  difference  in 
using  Dalton’s  law  of  mixtures  (for  the  mass  fraction  case)  instead  of  Amagat’s  law  (for 
the  volume  fraction  case).  Both  are  equivalent  and  both  give  the  same  answer.  For 
computational  purposes,  there  may  however  be  a  preference  for  one  over  the  other. 

The  third  way  to  express  the  mixture  properties  is  in  terms  of  the  mole  fraction, 

v  _  moles  of  vapor  _  molv 

A  y  "  -  ' 

mole  of  mixture  rnolmix 


The  mixture  density  is  then  obtained  as, 

_  M_mix¥v  _  kgy  moly  k^mix  _  molv 
Mv  kg  mix  kgv  mo^mix  mo^mL 
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where  Mmix  represents  the  mean  molecular  weight  of  the  mixture  and  Mv  is  the 
molecular  weight  of  the  vapor.  From  this  relation,  we  find, 


P  Yy 


=  avpv  = 


P  MyXy 

Mmix 


but 


P^v  _  kg  mix  k8v  m<->hix  _  hv  mo^mix  _  hv 

M mix  m3mix  molv  kgmix  n?mix  molv  ^ 


where  we  have  related  the  molecular  weight  to  the  volume.  Consequently,  the  mole 
fraction  is  analogous  to  the  volume  fraction. 

For  two-phase  flows  the  mass  fraction  and  the  volume  fraction  are  equivalent,  but 
generally  the  volume  fraction  works  better. 


P  -^vPv  O’  *V )p /  — 


hv  .  hi  _  hv  , 

3  3  +  3  3  3  + 

mrnix  mmix  ml  mmix 


k8l  _  kg  mix 
3  3 

mmix  mmix 


Fonmdation  for  Time-Accurate  Confutations:  The  above  development  applies  to  steady 
state  flows  only.  For  the  cavitation  problem,  it  is  imperative  that  we  consider  time- 
accurate  flows  to  be  able  to  simulate  the  unsteadiness  that  is  characteristic  of  flows  with 
attached  sheet  cavitation.  For  implementation  in  a  time-accurate  sense,  we  write  the 
dual-time  version  of  the  two-phase  equations.  We  include  the  physical  time  as  part  of  a 
four-dimensional  divergence  ( x,y,z,t )  operator  and  add  a  pseudo-time  term,  t,  for  the 
iterations  at  each  time  step  (or  for  iteration  to  the  steady  state  if  the  physical  time 
derivative  term  goes  to  zero).  We  then  use  conservation  variables, 

Q  =  (p,pw,pv.pw,p/i°  -  pj  in  the  physical  time  derivative  and  the  primitive  variables  in 
the  pseudo-time  derivative  along  with  a  preconditioning  matrix.  The  resulting  vector 
form  then  becomes: 


dQp  dQ  dE  dF  dG  rT 

—f  +  ^+— +  — +  — =  l7  +  v.r. 
or  at  ox  oy  oz 


where  H  represents  the  source  (and  sink)  terms  in  the  continuity  equations  representing 
phase  changes. 

The  two  continuity  equations  can  be  expressed  as  one  equation  for  the  liquid  and  one  for 
vapor  phase  or  as  the  global  continuity  equation  (the  sum  of  these  two)  and  either  of  the 
phasic  equations.  Here  we  express  the  system  in  terms  of  the  total  conservation  of  mass 
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equation  (the  sum  of  the  liquid  and  gas  phase  equations)  and  the  liquid  phase 
conservation  equation.  With  this  choice,  the  terms  in  the  equations  become, 
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The  system  is  completed  by  specifying  two  equations  of  state: 

Pl=Pl{p>T),  and  pv  =p v(p,T), 

given  earlier,  along  with  the  definition  of  the  mixture  density, 
p  =  azpz+(l-cxz)pv 

The  system  is  then  closed  by  specifying  analogous  thermodynamic  relations  for  the 
enthalpy  of  each  phase, 

h  =  hi  ( p,T )  and  hv=hv(p,T ) 

The  mixture  enthalpy  relation  likewise  follows  the  relation  for  the  total  density, 
ph  =  alplhl+(l-<xl)pvhv 


The  important  properties  that  appear  in  the  Jacobian  of  the  flux  term  and  also  in  the 
preconditioning  matrix  are  given  by  differentiating  the  density  function(s)  in  the 
continuity  equation(s)  to  get  the  partial  derivatives  of  density  with  respect  to  pressure, 
temperature  and  volume  fraction.  These  are  most  readily  obtained  by  computing  the 
Jacobian, 
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where  subscripts  refer  to  partial  differentiation.  For  example,  the  terms  with  the 
subscript,  p,  represent  derivatives  with  respect  to  the  pressure,  and  are  defined  as: 


Pip  ~ 


M 

dp 


Ju,T,  a 


Pvp  ~ 


(dPv 

<  dp 


)u,T,  a 


and 


Pv s 


dp 

V  dP  Ju,T,a 


\ 


=  «z 


dpi 

dp 


+(i 

Ju,T,a  1  & 


Ju,T,  a 


=  alPlp  +  (l  “  az  )p 


vp 


The  terms  with  the  subscript,  T,  represent  derivatives  with  respect  to  the  temperature. 
They  are  defined  as, 


Pit  = 


d_Pi 

l  dp  Jptu,a 


PvT  = 


^Pv"| 
dT  )piU,a 


and 


dT 


P  T  ~ 


)  p,u,a 


~alPlT  +(l~al)PvT 


p,u,  a 


Finally  the  derivative  with  respect  to  the  volume  fraction,  oq ,  also  appears.  This  term  is 


Pa  1 


dp_ 

da, 


pm,t 


=  P l  “Pv 


Note  that  derivatives  of  the  individual  densities  with  respect  to  the  volume  fraction  are 
identically  zero  since  the  phasic  densities  do  not  depend  upon  the  volume  fraction, 


fia 

da. 


\ 

)p,u,T 


=  0 


'*0 


da, 


Jp,u,T 


The  viscous  terms  in  the  two-phase  equation  system  are  identical  to  those  given  above 
except  for  the  additional  equation  and  are  not  re-written. 

The  primary  new  term  in  the  equations  for  the  two-phase  system  is  the  source  term,  H, 
which  ensures  local  conservation  of  mass  when  liquid  changes  to  vapor  or  when  vapor 
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changes  to  liquid.  These  source  terms  correspond  to  the  phase  transition  model.  For  the 
liquid  volume  fraction  equation,  the  source  term  includes  one  contribution  that  describes 
the  rate  of  conversion  of  vapor  to  liquid  and  one  for  the  conversion  of  liquid  to  vapor. 
The  conversion  from  vapor  to  liquid  appears  as  a  source  term  in  the  liquid  volume 
fraction  equation,  while  the  conversion  from  liquid  to  vapor  appears  as  a  source  term. 
These  terms  appear  on  the  right  hand  side  of  the  conservation  equation  and  take  the  form, 

„  da*  ,  dp/a z  ,  3pzazw  _  az  ,  (l-az) 

PIP  dx  dt  dx  ~  T i  Tv 


where  tz  and  xv  are  the  characteristic  times  for  conversion  from  liquid  to  vapor  and 

from  vapor  to  liquid  respectively.  These  times  are  chosen  as  being  proportional  to  the 
local  pressure  difference, 


for  p<pv: 

for  p>pv: 


l-o. 

l  _  i 

P~Pv 

H  1 

«-»l 

1 

c 

ty  T, 

q 

1  _  1 

P~Pv 

J_=o 

TZ  Xr 

q 

9 

and 


A  representative  choice  with  K  =  (Kv  / Kl)  =  l/(a[  +pv/pz)  is  shown  on  Fig.  2  for 

various  values  of  K  ranging  from  1.0  to  .001.  The  corresponding  ratio  of  the  speed  of 
sound  in  the  two-phase  mixture  to  the  speed  of  sound  in  the  liquid  is  given  on  Fig.  3.  This 
figure  shows  the  familiar  minimum  in  the  sound  speed  when  changing  from  100%  vapor 
to  100%  liquid.  Note  that  in  the  case  where  the  liquid  and  vapor  densities  are  equal,  that 
no  minimum  speed  of  sound  appears  between  the  liquid  and  gas  phases. 


Franc -Michel  Experiment:  Franc  and  Michel  [l,2]have  reported  a  series  of  cavitation 
experiments  on  hydrofoils  that  are  appropriate  for  calibration  and  validation  purposes  in 
the  above  models.  Their  experiments  show  that  turbulent  spots  from  exploding  nuclei 
can  remove  sheet  cavitation,  but  that  leading  edge  sheet  cavitation  is  resistant  to  the 
presence  of  free-stream  nuclei  and  to  the  boundary  layer  state,  whether  it  is  laminar  or 
turbulent.  Their  cavitation  results  show  that  the  cavitation  pattern  is  two  dimensional 
near  the  leading  edge  where  it  appears  at  large  angles  of  attack  and  high  cavitation 
numbers.  Mid-chord  cavitation  appears  at  middle  angles  of  attack  and  middle  values  of 
the  cavitation  number.  Mid-chord  cavitation  tends  to  be  three  dimensional.  When  the 
angle  of  attack  and  the  cavitation  number  are  both  low,  cavitation  occurs  at  the  trailing 
edge.  Trailing  edge  cavitation  is  generally  predominantly  two  dimensional. 

Careful  sequences  of  measurements  show  that  cavitation  moves  from  the  trailing  edge  to 
the  leading  edge  as  the  angle  of  attack  is  increased.  In  general,  attached  cavitation  is 
established  downstream  of  boundary  layer  separation,  but  as  the  number  of  free-stream 
nuclei  increases,  cavitation  moves  toward  the  minimum  pressure  point.  When  the 
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boundary  layer  transitions  from  the  laminar  to  the  turbulent  state,  the  turbulence 
generally  eliminates  attached  cavitation.  Finally,  Franc  and  Michel’s  experiments  show 
that  for  selected  angles  of  attack,  decreasing  the  cavitation  number  can  cause  cavitation  to 
first  appear,  then  disappear  and  finally  reappear  again.  This  non-monotonic  behavior  and 
particularly  the  fact  that  decreasing  the  cavitation  number  causes  cavitation  to  disappear 
suggests  that  non-equilibrium  effects  are  present. 


Predictions  from  Two-Phase  Model:  A  series  of  computations  has  been  computed  with 
the  dual-phase  model.  In  general,  these  results  are  unsteady  and  exhibit  fluctuations 
around  the  trailing  edge,  and  in  some  cases  the  leading  edge,  of  the  cavitation  region.  As 
the  cavitation  number  decreases  (and  the  size  of  the  cavitation  bubble  increases),  the 
flowfield  becomes  more  and  more  unsteady.  Some  representative  results  from  the  two- 
phase  cavitation  model  are  discussed  in  the  present  section.  Most  of  the  results  shown 
have  been  computed  on  a  multi-block  grid  involving  typically  20  to  30  blocks,  although 
where  noted  the  results  have  been  obtained  on  a  single-block  grid.  Representative  grids 
are  shown  on  Figs.  6  and  7  for  both  the  near  field  and  the  far  field.  The  plots  that  show 
the  multi-block  results  include  the  outline  of  the  multiple  blocks.  The  present  two  figures 
are  for  grids  of  30,000  nodes  each.  The  multi-block  configuration  in  Fig.  6  uses  30 
blocks. 

As  a  first  example  of  the  predictions  from  the  two-phase  model,  we  show  results  for  flow 
over  the  NACA  66  (mod)  hydrofoil  used  in  the  Shen-Demotakis  experiment  referred  to 
above.  Figure  8  shows  the  density  and  pressure  contours  for  flow  at  4  degrees  angle  of 
attack  and  a  cavitation  number  of  0.91.  The  Reynolds  number  is  one  million,  and  the 
density  ratio  between  liquid  and  vapor  was  taken  as  100.  Since  the  density  is  constant 
over  the  entire  liquid  region,  the  density  contours  appear  only  in  the  cavitation  bubble 
region.  The  pressure  contours,  however,  are  distributed  over  the  entire  flowfield  and,  in  a 
qualitative  sense,  retain  the  familiar  pressure  distribution  of  single-phase  flows.  Close 
inspection  of  the  flow  near  the  cavitation  region,  however,  shows  the  pressure  is 
essentially  constant  inside  the  cavity  as  is  frequently  assumed  in  simple  cavitation 
models.  We  note  that  the  pressure  in  this  region  does  retain  weak  variations  rather  than 
being  strictly  constant.  In  addition,  the  (nearly)  constant  pressure  region  appears  as  a 
result  of  the  computation,  not  as  the  result  of  an  assumption  in  the  modeling.  Again,  the 
background  segments  in  the  plot  represent  boundaries  of  the  multi-blocks,  and  are  not  a 
part  of  the  solution. 

A  close-up  view  of  the  details  in  the  cavitation  region  for  this  same  computation  is  shown 
in  Fig.  9.  The  velocity  vectors  in  this  region  show  the  presence  of  a  re-entrant  jet  on  the 
downstream  end  of  the  cavitation  region.  This  narrow  re-entrant  jet  flows  forward 
against  the  free-stream  flow  direction  and  lies  between  the  rear-most  part  of  the  cavity 
and  the  surface  of  the  hydrofoil.  In  addition,  the  flowfield  gives  the  suggestion  that  a 
portion  of  the  cavity  is  about  to  be  shed  into  the  wake.  Long-term  computations  do  show 
that  the  cavity  length  fluctuates  slowly  as  the  result  of  such  periodic  shedding. 
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Figure  10  shows  the  pressure  contours  over  a  wedge  patterned  after  experiments  by 
Ceccio.  This  figure  shows  the  pressure  contours  on  the  fully  wetted  wedge  as  a  function 
of  the  distance  along  the  wedge  surface.  The  results  contrast  the  pressure  contours  when 
cavitation  is  present  with  those  when  the  wedge  surface  is  fully  wetted.  The  presence  of 
cavitation  shows  a  decided  change  in  the  surface  pressure  that  should  be  measurable 
experimentally. 

The  next  several  figures  concern  the  flow  over  a  NACA  16-012  hydrofoil  like  that  used 
in  the  Franc  and  Michel  experiments.  Both  multi-block  and  single-block  grid 
computations  were  obtained  for  this  airfoil.  In  both  cases,  the  grid  contained  from  20,000 
to  50,000  grid  points.  The  intent  of  these  computations  was  to  provide  a  fine  enough  grid 
structure  to  ensure  that  the  conclusions  deduced  from  the  results  were  not  grid  dependent. 
The  following  figures  show  representative  solutions  on  these  and  similar  grids. 

Fig.  11  shows  an  attempt  at  a  laminar  flow  solution  over  the  NACA  16-012  hydrofoil. 
The  flow  conditions  of  interest  were,  of  course,  turbulent,  but  a  series  of  laminar  flow 
solutions  were  attempted  to  understand  the  underlying  flowfield  without  the  uncertainties 
introduced  by  a  turbulence  model.  For  all  but  the  very  lowest  Reynolds  number 
conditions  attempted,  the  laminar  solutions  resulted  in  laminar  separation  at  the  trailing 
edge,  and  the  resulting  flowfield  was  highly  unsteady.  Figure  1 1  shows  the  flowfield  at 
one  instant  of  time  in  a  computation  at  a  Reynolds  number  of  75,000.  The  highly 
unsteady  flowfield  is  clearly  indicated  here. 

The  emphasis  in  modeling  cavitation  is  on  finding  hydrofoils  that  contain  laminar  flow 
over  the  leading  edge  where  the  cavity  develops,  but  turbulent  flow  over  the  trailing  edge 
to  prevent  laminar  separation  and  to  stabilize  the  flowfield.  The  present  laminar 
computations  demonstrated  amply  that  the  unsteadiness  from  a  completely  laminar 
solution  would  not  allow  cavitation  studies  to  be  done  without  a  turbulence  model.  It  is 
imperative  in  cavitation  studies,  however,  that  the  turbulence  model  provide  a  non- 
cavitating  flowfield  that  is  qualitatively  proper.  Specifically  this  means  that  the  flow 
must  be  laminar  over  the  leading  edge  of  the  hydrofoil  to  allow  cavitation,  but  turbulent 
at  the  trailing  edge  to  prevent  separation.  Experimental  evidence  shows  that  the  presence 
of  turbulent  flow  near  the  leading  edge  prevents  the  formation  of  attached  sheet 
cavitation.  The  present  model  is  probably  not  sensitive  to  the  presence  of  turbulence  in 
establishing  cavitating  regions,  but  it  is  important  to  establish  a  correct  flowfield  to 
enable  later  improvements. 

One  of  the  characteristics  of  the  Franc-Michel  experiments  is  the  manner  in  which  the 
cavitation  position  moves  from  the  trailing  edge  to  the  leading  edge.  Franc  and  Michel 
also  noted  that  the  locations  of  transition  to  turbulence  and  separation  likewise  varied 
with  angle  of  attack.  The  present  results  are  computed  with  the  k-e  model,  and  the  free- 
stream  turbulence  has  been  carefully  adjusted  to  get  the  transition  location  and  the 
turbulence  behavior  to  match  that  of  the  experiments.  Results  for  two  different  levels  of 
free-stream  turbulence  are  illustrated  below  to  indicate  the  manner  in  which  the  flowfield 
changes  with  the  turbulence  level,  and  the  sensitivities  involved  in  matching  the  fully 
wetted  airfoil  experiments. 
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As  a  first  example  of  the  effects  of  a  turbulence  model,  we  present  on  Fig.  12  the  velocity 
contours  and  the  ratio  of  the  eddy  viscosity  to  the  laminar  viscosity  for  the  NACA  16-012 
hydrofoil  at  zero  degrees  angle  of  attack.  The  present  computation  uses  third-order 
upwind  differencing  and  implicit  discretization  in  time  and  is  for  a  cavitation  number  at 
which  the  flow  is  fully  wetted  and  cavitation  does  not  occur.  As  can  be  seen  from  the 
figure,  the  solution,  both  in  terms  of  the  velocity  and  the  turbulence  level  is  symmetric 
about  the  hydrofoil  as  would  be  expected.  The  level  of  free-stream  turbulence  in  this 
zero-degree  condition  has  been  adjusted  to  match  the  experimental  observations  at  this 
flow  condition.  Hereafter,  we  refer  to  this  free  stream  turbulence  level  as  the  'low' 
turbulence  level.  As  seen  in  the  figure,  the  effects  of  turbulence  are  significant  only  in 
the  boundary  layers  near  the  trailing  edge  (approximately  the  last  10%)  of  the  hydrofoil 
and  in  the  wake  behind  the  hydrofoil.  This  local  turbulent  region  is,  however,  sufficient 
to  remove  the  large-scale  unsteadiness  that  was  observed  in  the  laminar  calculation  at  Re 
=  75,000  (Fig.  10).  In  contrast  to  that  laminar  case,  the  present  turbulent  case  converged 
very  well  to  a  steady  state  solution  with  a  very  small  recirculation  zone  at  the  trailing 
edge.  Again,  both  the  steady  solution  and  the  presence  of  turbulence  over  the  aft  10%  of 
the  hydrofoil  are  in  agreement  with  the  Franc  and  Michel  measurements. 

Corresponding  results  for  the  NACA  16-012  hydrofoil  at  three  degrees  angle  of  attack  are 
shown  on  Figs.  13  and  14.  The  chord  Reynolds  number  for  this  case  is  again  300,000 
and  the  same  grid  and  solution  procedures  were  used.  Figure  13  shows  contours  of  the 
velocity  and  the  ratio  of  turbulent  to  laminar  viscosity  for  this  three-degree  angle  of 
attack.  Here  the  velocity  clearly  shows  non-symmetric  (upper-to-lower  surface)  contours 
because  of  the  angle  of  attack.  The  eddy  viscosity  likewise  shows  a  similar  asymmetry. 
The  interesting  fact  is  that  turbulence  starts  at  about  mid-chord  on  the  upper  surface, 
while  its  transition  region  on  the  lower  surface  is  closer  to  the  trailing  edge  than  in  the 
zero-degree  angle  of  attack  case. 

The  corresponding  pressure  contours  for  the  three-degree  angle  of  attack  case  are  given 
on  Fig.  13  along  with  a  plot  of  the  multi-block  grid  that  was  used  for  the  computations. 
The  computed  pressure  distribution  on  the  surface  is  in  good  agreement  with 
experimental  measurements  at  these  conditions.  Again,  the  flowfield  at  a  =  3  degrees 
results  in  a  well  converged,  steady  solution.  Finally,  a  near-field  view  of  the  ratio  of 
turbulent  to  laminar  viscosity  is  shown  on  Fig.  15.  More  details  of  the  asymmetry  in  this 
figure  are  discussed  below. 

The  effect  of  the  freestream  turbulence  level  on  the  flowfield  solution  is  shown  on  Figs. 
16  and  17  for  two  different  free-stream  turbulence  levels.  Figure  16  shows  the  solution 
for  the  ’low'  freestream  turbulence  conditions  of  the  results  in  Figs.  12  to  15.  The  top  half 
of  Fig.  16  shows  a  near-field  view  of  the  velocity  vectors  over  the  rear  portion  of  the 
airfoil  for  the  zero-degree  angle  of  attack,  while  the  bottom  half  shows  the  results  for  the 
three-degree  case.  Similar  results  are  shown  on  Fig.  17  for  a  high  free  stream  turbulence 
level. 
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Changing  the  free-stream  turbulence  level  for  the  zero-degree  angle  of  attack  case  (the 
top  plots  of  Figs.  16  and  17)  has  very  little  effect  on  the  results.  The  velocity  vectors  in 
both  cases  are  quite  similar.  (Note  that  the  A  different  (single  block)  grid  was  used  for 
the  high  turbulence  results  on  Fig.  17  (this  grid  is  given  in  Fig.  7),  so  that  the  locations  of 
the  velocity  vectors  is  somewhat  different,  but  overall  the  two  results  are  qualitatively 
similar. 

Results  for  the  three-degree  angle  of  attack  are,  however,  considerably  different.  The 
high  free-stream  turbulence  calculation  (Fig.  17)  shows  a  thick  turbulent  boundary  layer 
on  the  upper  side  of  the  airfoil  and  a  thin  one  on  the  lower  side,  while  the  low  free-stream 
turbulence  case  (discussed  above)  shows  a  thicker  boundary  layer  on  the  lower  surface. 
Close  inspection  of  this  result  shows  that  the  reason  that  the  boundary  layer  is  thicker  on 
the  lower  surface  is  that  there  is  a  separation  region  at  about  85%  chord.  Comparison 
with  the  turbulence  profiles  in  16  indicates  that  this  is  a  laminar  separation  and  that  the 
turbulence  appears  downstream  of  the  laminar  separation  position.  This  separated  region 
is  the  reason  that  the  lower  surface  has  a  thicker  boundary  layer.  The  low  turbulence 
results  are  in  agreement  with  the  observations  of  Franc  and  Michel  indicating  that  this 
prediction  is  qualitatively  in  agreement  with  the  experiment  while  the  results  at  the  higher 
turbulence  level  (which  is  more  representative  of  ‘typical’  conditions)  are  not. 

The  pressure  distributions  on  the  hydrofoil  for  these  two  cases  are  shown  on  Fig.  18.  The 
change  in  the  turbulence  level  causes  some  changes  in  the  pressure  distribution  near  the 
aft  end  of  the  hydrofoil  where  the  flow  separates  from  the  surface. 

Finally,  a  comparison  with  the  results  of  Franc  and  Michel  is  given  on  Fig.  19.  In  this 
figure,  the  locations  of  several  boundary  layer  events  are  tracked  as  a  function  of  the 
angle  of  attack  between  negative  ten  and  positive  ten  degrees.  In  particular,  the  plot 
shows  the  location  of  the  laminar  separation  point,  the  start  of  transition  to  turbulence,  the 
end  of  transition  and  the  location  of  turbulent  separation  on  the  upper  side  of  the  airfoil. 
As  an  example,  the  laminar  separation  point  on  the  upper  surface  lies  at  about  90%  chord 
when  the  angle  of  attack  is  minus  ten  degrees.  As  the  angle  of  attack  is  increased,  the 
laminar  separation  point  moves  forward  in  an  approximately  linear  fashion  until  the  zero- 
degree  angle  is  reached.  At  zero  degrees,  the  location  of  laminar  separation  is  at  about 
80%  chord.  As  the  angle  of  attack  becomes  positive,  the  point  of  laminar  separation 
continues  to  move  forward,  but  is  very  rapidly  overtaken  by  the  start  of  transition  to 
turbulence  at  which  point  laminar  separation  ceases  to  exist. 

Representative  predictions  from  the  present  computations  are  included  on  Fig.  19  at  -3,  0, 
and  +3  degrees  angle  of  attack.  The  location  of  the  laminar  separation  point  is  seen  to 
track  the  France  and  Michel  results  quite  well  at  the  -3  and  the  0  degree  locations.  At  the 
+3  angle  of  attack,  laminar  separation  is  not  observed  on  the  upper  surface,  again  in 
agreement  with  the  experimental  results. 

Results  for  the  'start'  of  transition  are  also  indicated  for  the  three  angles  of  attack.  The 
computational  results  for  the  start  of  transition  are  seen  to  track  the  Franc/Michel  results 
qualitatively,  although  they  are  some  5%  downstream  of  the  Franc/Michel  results.  This 


difference  in  location  probably  arises  because  of  different  methods  for  identifying  the 
start  of  transition.  The  k-e  model  does  not  provide  a  direct  measure  of  'transition'  and  it 
must  be  inferred  from  the  velocity  profiles.  The  computed  results  do  show  that  the 
transition  location  moves  forward  gradually  in  the  negative  angle  of  attack  region,  and 
then  moves  much  more  rapidly  in  the  positive  angle  of  attack  region. 

Finally,  the  computations  show  a  single  turbulent  separation  point  at  the  positive  angle  of 
attack  location.  This  turbulent  separation  point  occurs  at  about  95%  chord.  Again,  the 
absence  of  turbulent  separation  (on  the  upper  side)  at  the  -3  and  0  degree  angle  of  attack 
cases  is  in  agreement  with  the  experimental  observations.  Overall,  the  present  results 
indicate  that  the  turbulence  model  can  reproduce  the  observed  turbulence  conditions  on 
the  hydrofoil  in  the  fully  wetted  condition  and  provide  an  appropriate  background  for 
employing  a  cavitation  model. 

A  substantial  amount  of  cavitation  results  from  the  single-phase  and  the  two-phase 
models  has  been  presented  in  Ref.  3  and  shows  that  the  two  models  give  reasonable 
predictions  of  cavitation  for  various  conditions.  The  two-phase  model  provides 
capability  for  the  quasi-equilibrium  results  observed  in  numerous  experiments  and  so  is 
the  preferred  model  for  additional  calculations. 
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Fig.  1 .  Density-pressure  curves  in  transition  region  for  continuous  equation  of  state  used 
in  the  single-fluid  model.  A p  =  1 00,  200,  300  and  400  N/m2. 
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Fig.  2.  Logarithm  of  the  speed  of  sound  of  continuous  equation  of  state  as  a  function  of 
pressure  in  transition  regime.  A p  =  100,  200,  300  and  400  N/m2. 
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Fig.  3.  Mach  number  contours  for  flow  over  an  ellipse  at  -2 

degrees  angle  of  attack  based  on  single-phase  cavitation 
model.  Transonic  flow  region  outlines  cavitation  zone. 


Fig.  4.  Ratio  of  Kv  /  Kt  for  various  vapor  to  liquid  density  ratios  as  a 
function  of  volume  fraction  of  liquid. 


Speed  of  Sound  in  Mixture 


o  of  speed  of  sound  in  two-phase  mixture  to  sound  speed  in 
id  as  a  function  of  volume  fraction  of  liquid  for  various  ratios 
apor-to-liquid  density  ratios. 
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Fig.  8.  Computed  density  and  pressure  contours  on  a  NACA  66  hydrofoil  based  upon 
the  two-phase  model,  a  =  4  degrees;  a  =  0.9 1 ;  Re  =  1*1 06;  pL  /p  =  1 00. 
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Fig  9.  Near-field  view  of  cavitation  region  in  Fig.  5  including  velocity  vectors  and 
density  contours. 


Fig.  10.  Computed  cavitating  and  non-cavitating  wall  pressure  distribution 
26.6  degree  wedge. 


Fig.  11.  Instantaneous  snapshot  of  stream-wise  velocity  contours  for  laminar  flow 
past  a  NACA  16-012  hydrofoil  at  a  chord  Reynolds  number  of  75,000. 


Fig.  12.  Velocity  contours  and  the  ratio  of  eddy  viscosity  to  laminar  viscosity  for  the 
NACA  16-012  hydrofoil.  Zero  degrees  angle  of  attack;  Re  =  300,000. 


NACA  16-012  hydrofoil.  Three  degrees  angle  of  attack;  Re  =  300,000. 
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3rd-order  upwind  biased  differencing,  2-eq  Turijulence  model 
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Pressure  Coefficients  for  N  ACA 1 6-01 2 


Fig.  18.  Single  block  grid  used  for  computation  of  high  turbulence  level  solution. 


Attached  cavitation  and  the  boundary  layer 


Position  of  boundary  layer  characteristics  versus  angle  of  attack. 
NACA  16-012  hydrofoil;  Upper  side. 


